
********************************************************************************
** This code replicates our pooled analysis using the data used in Part 2 of
** Levinson (AER, 2016). The results are displayed in the first 7 columns of
** Table A15 in our Appendix as well as in Figure 7.

			 
#delimit ;
clear all;
set more 1;


*---------------------------------------------------;
*  Get the data created by Part2data.do;
*---------------------------------------------------;

	use "Data/Part2.dta", clear;
	
	merge m:1 combinedID 
		using "Data/Part1_NoEnergy_1.dta", generate(part2merge);
	
		gen kBTUmonthelec = 3412.14163 * kwhmonth /1000; *Note, KBTU not MBTU;
		gen lnmonthelec = 100*ln(kwhmonth); *Note BTU not KBTU or MBTU;
		gen kBTUmonthgas = 99976.129 * ngmonth / 1000 ;
		gen lnBTUmonthgas = ln(kBTUmonthgas*1000); *Note BTU not KBTU or MBTU;
		
		*weather data is in 10ths of CDDs.  Convert to  degrees;
		  foreach VARIABLE of varlist cddzip20 cddzip30 cddzip40 hddzip20 hddzip30 hddzip40 avgCDDzip20 avgCDDzip30 avgCDDzip40 avgHDDzip20 avgHDDzip30 avgHDDzip40 {;
		   	  replace `VARIABLE' = (9/5)*`VARIABLE' /10;
		  };
		
		keep if dwltype==1 & constYRS ~=97 & elecwater==0 & elecheat==0;
		
		*drop utility outliers. Chong drops <2kWh/day or >80kWh/day  (4% of sample);
		*I include more;
		
			gen over = (kwhmonth>3000);
			gen under = (kwhmonth<100);
			
			summ over under;
			drop if over | under;
		
		drop if constYRS==97;
		
		gen lnrescnt = ln(rescnt);		
		gen lnsqftK = ln(sqftK);
		
		gen roomair = (actyp1==1) | (actyp2==1) | (actyp3==1);
		gen centralair = (ctlacage>=1) & (ctlacage<=3);
		gen sqftKxcentralair = sqftK*centralair;
		gen sqftKxelecheat = sqftK*elecheat;

		gen lnsqftKxcentralair = lnsqftK * centralair;
		gen lnbuildingage = ln(buildingage);
		
		gen lnyrs_res = ln(yrs_res);
		
		drop if hohblk1==97; *head of household black missing;
		
		gen year09 = (year>=2005);
		gen remodeled = (remod==2);
		
		
		preserve;
			collapse czt24, by(servzip);
			rename czt24 climatezone;
			save "Data\Temp\cz", replace;
		restore;
		merge m:1 servzip using "Data\Temp\cz";
		drop _merge;
		
		
	save "Data\Temp\tempandenergy", replace;

*---------------------------------------------------;
*    Predict elect as a fn of cdd and hdd;
*---------------------------------------------------; 

	use "Data\Temp\tempandenergy", clear;
	
		*Express degree days in 10s;
		
		local UNITS 1;
		
		gen elecboth = 0;  replace elecboth = 1 if elecstove==1 & elecoven==1;
		
		egen houseFE = group(combinedID);
			
		*replace degree days with 10s of degree days;
		*weather data is converted to degrees, above.  Convert to 10s of degrees for presentation of regression coefficients;;
		  foreach VARIABLE of varlist cdd65 cddzip20 cddzip30 cddzip40 hdd65 hddzip20 hddzip30 hddzip40 avgCDDzip20 avgCDDzip30 avgCDDzip40 avgHDDzip30 {;
		   	  gen `VARIABLE'_0 = `VARIABLE' /`UNITS';
		  };
	
		*Interact air conditioning with CDD;
			gen ACxCDD = ((centralair==1)|(roomair==1)) * cddzip30;
			


	*generate cdd interaction and hdd interaction ;
	*omit 1978-82, when constYRS=7);
		drop if constYRS==97;
		
		foreach VIN of numlist  1 2 3 4 5 7 8 9 10 11 12 {;
			gen cddXvintage`VIN' = cddzip30_0*(constYRS==`VIN');
			gen hddXvintage`VIN' = hddzip30_0*(constYRS==`VIN');
						
			gen cddXvinXnoAC`VIN' = cddXvintage`VIN' * max(0,1-centralair-roomair);  *interact with "no AC" so regular interaction coefficients are the houses with AC;
		};	


		forvalues VIN = 1/38 {;
			gen cddXbuiltyr`VIN' = cddzip30_0*(builtyr==`VIN') if year09==0;
			gen hddXbuiltyr`VIN' = hddzip30_0*(builtyr==`VIN') if year09==0;
		};	

	
		*Add age of heater and AC, and interact cdd and hdd with those;	
	
			foreach AGE of numlist   1 2 3 4 5   {;
				gen cddXacAGE`AGE' = cddzip30_0*(clcntage==`AGE');
				gen hddXhhAGE`AGE' = hddzip30_0*(htsysage==`AGE');
			};	
	
	egen countyFE = group(county);
		
	*generate climatezone-month fe's;
		egen climatezonemonthFE = group(cecfast month);
	*generate county-month fe's;
		egen countymonthFE = group(county month);	
	*generate zip-month fe's;
		egen zipmonthFE = group(servzip month);


	save "Data\Temp\Part2regs", replace;

	*--------------;
	*ELECTRICITY;
	*--------------;
	
	local CONSTRUCTFE const4049 const5059 const6069 const7074 
		const7577 const7882 const8392 const9397 const9800 const0104 const0508;
	local housechars  lnsqftK numroom elecboth remodeled;
	local reschars lnyrs_res lnrescnt lnhhinc nr0_5  nr65_99  college anydisabled hohblk1 hohlat1 own	;
	local appliances rfnum roomair centralair ; *drop freezers (never significant);
	
	local FUEL 	lnmonthelec;
	local ALTERNATIVEFUEL kwhmonth;
	local DEC 3;
	local OUTFILE ""Tables\Replications.doc"";

	local CLUSTER countyFE;

	use "Data\Temp\Part2regs", clear;

	
	*First get number of houses represented, by year;
		quietly reg `FUEL' cddzip30_0 `housechars' `reschars' `appliances' cddXvintage* year09 `CONSTRUCTFE' countymonthFE countyFE;
		gen sample=e(sample);
		
		egen HnumELEC2002 =group(combinedID ) if year09==0 & sample;
		egen HnumELEC2009 =group(combinedID ) if year09==1 & sample;

		summ HnumELEC*;

	gen summ=month>4&month<10&(year==2002|year==2008);
	gen aa=cddzip30_0*summ;
	egen annCDD1=sum(aa), by(combinedID);
	egen annCDD=mean(annCDD1), by(servzip);
	drop summ aa annCDD1;
	
	gen CDDxsqft=cddzip30_0*(sqftK-1.75);
	
	gen annCDDxCDD1=cddzip30_0*(annCDD<=400);
	gen annCDDxCDD2=cddzip30_0*(annCDD>400)*(annCDD<=800);
	gen annCDDxCDD3=cddzip30_0*(annCDD>800)*(annCDD<=1200);
	gen annCDDxCDD4=cddzip30_0*(annCDD>1200);
	
	encode combinedID, gen(premiseID);
	xtset premiseID yearmo;

	foreach cz in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16  {;
		gen hdd_`cz' = hddzip30_0*(climatezone==`cz');
		gen cdd_`cz' = cddzip30_0*(climatezone==`cz');
		gen CDDxsqftx_`cz'=CDDxsqft*(climatezone==`cz');
	};

	gen cddXbuiltyr7879=cddXbuiltyr13+cddXbuiltyr14;
	gen cddXbuiltyr8082=cddXbuiltyr15+cddXbuiltyr16+cddXbuiltyr17;
	
	rename cddXbuiltyr12 cddXbuilt_12;
	//drop if climatezone==.;
		
********************************************************************************
** The following models produce the results displayed in the first 7 columns of
** Table A15. The final model (model G) produces the results that are displayed
** in Figure 7 -- the results corresponding to Levinson's data.

	*column A -- Log, All Years;
	local X A;
	areg lnmonthelec cddzip30_0 avgCDDzip30_0 `CONSTRUCTFE'  `housechars' `reschars' `appliances' cddXvintage* year09 i.month , absorb(countyFE) cluster(`CLUSTER') ;
		outreg2 cddzip30_0 avgCDDzip30_0 hddzip30_0  `CONSTRUCTFE'  `housechars' `reschars' `appliances' cddXvintage* year09 i.month using 
			`OUTFILE',  slow(1000) se alpha(.05) symbol(*) replace dec(3) ctitle(log) ;

	*column B -- Levels, All Years;
	local X B;
	areg kwhmonth cddzip30_0 avgCDDzip30_0 `CONSTRUCTFE'  `housechars' `reschars' `appliances' cddXvintage* year09 i.month , absorb(countyFE) cluster(`CLUSTER') ;
		outreg2 cddzip30_0 avgCDDzip30_0 hddzip30_0  `CONSTRUCTFE'  `housechars' `reschars' `appliances' cddXvintage* year09 i.month using 
			`OUTFILE',  slow(1000) se alpha(.05) symbol(*) append dec(`DEC') ctitle(level) ;

	*column C -- Levels, noavgCDD;
	local X C;
	areg kwhmonth cddzip30_0  hddzip30_0 `CONSTRUCTFE'  `housechars' `reschars' `appliances' cddXvintage* year09 i.month , absorb(countyFE) cluster(`CLUSTER') ;
		outreg2 cddzip30_0 avgCDDzip30_0 hddzip30_0  `CONSTRUCTFE'  `housechars' `reschars' `appliances' cddXvintage* year09 i.month using 
			`OUTFILE',  slow(1000) se alpha(.05) symbol(*) append dec(`DEC') ctitle(noavgCDD) ;

	*column D -- Levels, FE;
	local X D;
	areg kwhmonth cddzip30_0  hddzip30_0 `CONSTRUCTFE'  `housechars' `reschars' `appliances' cddXvintage* year09 i.month , absorb(premiseID) cluster(`CLUSTER') ;
		outreg2 cddzip30_0 avgCDDzip30_0 hddzip30_0  `CONSTRUCTFE'  `housechars' `reschars' `appliances' cddXvintage* year09 i.month using 
			`OUTFILE',  slow(1000) se alpha(.05) symbol(*) append dec(`DEC') ctitle(FE) ;

	*column E -- Levels, 78-82;
	local X E;
	areg kwhmonth cddzip30_0  hddzip30_0  `CONSTRUCTFE'  `housechars' `reschars' `appliances' cddXvintage* year09 i.month if (constYRS==6|constYRS==7), absorb(premiseID) cluster(`CLUSTER') ;
		outreg2 cddzip30_0 avgCDDzip30_0 hddzip30_0  `CONSTRUCTFE'  `housechars' `reschars' `appliances' cddXvintage* year09 i.month using 
			`OUTFILE',  slow(1000) se alpha(.05) symbol(*) append dec(`DEC') ctitle(75-82) ;

	*column F -- Levels, 2003;
	local X F;
	areg kwhmonth cddzip30_0  hddzip30_0  `CONSTRUCTFE'  `housechars' `reschars' `appliances' cddXvintage* year09 i.month if (constYRS==6|constYRS==7)&year09==0, absorb(premiseID) cluster(`CLUSTER') ;
		outreg2 cddzip30_0 avgCDDzip30_0 hddzip30_0  `CONSTRUCTFE'  `housechars' `reschars' `appliances' cddXvintage* year09 i.month using 
			`OUTFILE',  slow(1000) se alpha(.05) symbol(*) append dec(`DEC') ctitle(2003) ;

	*column G -- Levels, Years;
	local X H;
	areg kwhmonth cddzip30_0  hddzip30_0  `CONSTRUCTFE'  `housechars' `reschars' `appliances' cddXbuiltyr* year09 i.month if (constYRS==6|constYRS==7)&year09==0, absorb(premiseID) cluster(`CLUSTER') ;
		outreg2 cddzip30_0 avgCDDzip30_0 hddzip30_0  `CONSTRUCTFE'  `housechars' `reschars' `appliances' cddXvintage* year09 i.month using 
			`OUTFILE',  slow(1000) se alpha(.05) symbol(*) append dec(`DEC') ctitle(years) ;
